rm(list = ls(all = TRUE))

gc()
gc()

seed = 123456
set.seed(seed)

1 libraries

library(scales)
library(ggplot2)
# packageVersion("ggplot2")
library(magrittr)
# library(crayon)
library(rstatix)
# library(ggpubr)
# library(Compositional)
# library(equalCovs)
# library(psych)
library(corrplot) 
# library(pdist)
library(RColorBrewer)
library(plotly)
library(heatmaply)
# library(Hmisc)
library(htmlwidgets)
library(viridis)
library(ggh4x)


`%nin%` = Negate(`%in%`)

2 f-ctions

3 predefined contrasts

fpi = file.path('..', 'input')
fn = 'contrasts.txt'

comp = data.table::fread(file.path(fpi, fn), header = FALSE)

4 metabolomics

dir.create(file.path('..', 'reports', 'metabolomics'))

cnt = 1
fpi = file.path('..', 'input')


fn = "data_targeted.xlsx"
omics = openxlsx::read.xlsx(file.path(fpi, fn), sheet = 'metabolomics-leaves')
data.table::setDF(omics)

Genotype = 'Desiree'
omics = cbind(Genotype, omics)


PostInductionDay = as.numeric(stringr::str_split_i(gsub('S', '', omics$SampleID), '_', i = 2))
omics = cbind(PostInductionDay, omics)

omics = omics[omics$PostInductionDay %in% 1:14, ]

Condition = gsub('_.*', '', omics$SampleID)
omics = cbind(Condition, omics)
omics = omics[-which(omics$Condition == 'W' & omics$PostInductionDay %in% c(8,14)), ]

omics$Condition = factor(omics$Condition, levels = c('C', 'H', 'D', 'HD', 'W'))

table(omics$Condition, omics$PostInductionDay)
##     
##      1 7 8 14
##   C  4 4 4  4
##   H  4 4 4  4
##   D  0 0 4  4
##   HD 0 0 4  4
##   W  4 4 0  0
mygroup = paste(omics$Condition, stringr::str_pad(omics$PostInductionDay, 1, pad = "0"), sep = '_')
# table(mygroup %in% comp$V1)
# table(mygroup %in% comp$V3)

omics = cbind(mygroup, omics)

4.1 palette

n = 6


# my.ggplot.palette = brewer.pal(8, "Set2")[c(8, 2, 6, 4, 5, 3)]#[c(8, 2, 6, 7, 5, 3)]
my.ggplot.palette = c('#7F7F7F', '#C00000', '#FFC000', '#ED7D31',  
                      # '#70AD47', 
                      '#4472C4')
pie(rep(1, length(my.ggplot.palette)), 
    col = my.ggplot.palette, 
    main = 'ggplot palette',
    labels = c('C', 'H', 'D', 'HD', 
               # 'HDW', 
               'W'))

extend.palette = cbind(my.ggplot.palette, c('C', 'H', 'D', 'HD', 
                                            # 'HDW', 
                                            'W'))
colnames(extend.palette) = c('colour', 'group')
extend.palette = merge(extend.palette, omics[, 1:4], by.x = 'group', by.y = 'Condition', all = TRUE)


par(mfrow = c(2,2))

nn = 9
# (intervals = cut(c(-1,1), n*2))
intervals = cut(c(-1,1), 
                breaks = setdiff(seq(-1, 1, 0.1), 0), 
                include.lowest = TRUE)

my.heatmap.col1 = c(rev(brewer.pal(nn, 'Reds')), 'white', brewer.pal(nn, 'Blues'))
# my.heatmap.col1 = cividis(19, direction = -1)
pie(rep(1, length(my.heatmap.col1)), 
    col = rev(my.heatmap.col1), 
    labels = levels(intervals), 
    main = 'heatmap palette 1')


my.heatmap.col2 = c(my.heatmap.col1[1:3], 
                    rep('grey90', length(4:(nn*2-2))),
                    my.heatmap.col1[(nn*2-1): (nn*2 + 1)])
pie(rep(1, length(my.heatmap.col1)), 
    col = rev(my.heatmap.col2), 
    labels = levels(intervals), 
    main = 'heatmap palette 2')


par(mfrow = c(1,1))

4.2 density

# column in which measurements start
nx = 6

data = omics
# colnames(data)[nx:ncol(data)]

data[, nx:ncol(data)] = apply(data[, nx:ncol(data)] , 2, as.numeric)

# https://www.cookbook-r.com/Manipulating_data/Converting_data_between_wide_and_long_format/
data.long = tidyr::gather(data, molecule, measurement, colnames(data)[nx]:colnames(data)[ncol(data)], factor_key=FALSE)


data.long$measurement = as.numeric(data.long$measurement)

data.long$molecule = factor(data.long$molecule, levels = colnames(data)[nx:ncol(data)])

data.long$PostInductionDay = ordered(data.long$PostInductionDay, levels = sort(unique(data.long$PostInductionDay )))



data.long = data.long[!is.na(data.long$measurement), ]


# X11()
ggplot(data.long[data.long$Genotype == 'Desiree', ], aes(log2(measurement), fill = Condition, colour = Condition)) +
  geom_density(alpha = 0.1) + 
  ggtitle('Desiree')  +
  theme_classic() +
  scale_colour_manual(name = 'Legend',
                      values = my.ggplot.palette) +
  scale_fill_manual(name = 'Legend',
                    values = my.ggplot.palette) 

# X11()
ggplot(data.long[data.long$Genotype == 'Desiree', ], aes(log2(measurement), fill = mygroup , colour = mygroup )) +
  geom_density(alpha = 0.1) + 
  ggtitle('Desiree')  +
  theme_classic() +
  scale_colour_manual(name = 'Legend',
                      values = unique(extend.palette[, 2:3])[, 1]) +
  scale_fill_manual(name = 'Legend',
                    values = unique(extend.palette[, 2:3])[, 1]
                    ) 

data.long.SE = summarySE(data = data.long, 
                         measurevar="measurement", 
                         groupvars=c("molecule", "Genotype", "Condition", "PostInductionDay", "mygroup"), 
                         na.rm = TRUE)
ncol = 6

jitter = position_jitter(width = 0.0, height = 0.0)
dodge = position_dodge(width=0.0)

4.3 plot measurements

cd = sort(unique(data.long.SE$PostInductionDay))

mypallette = my.ggplot.palette

# useful in case of more genotypes
lims = data.long.SE %>% 
  dplyr::group_by(molecule) %>% 
  dplyr::summarise(Panel = molecule, 
                   ymin=min(mean-se),ymax=max(mean+se),
                   n = 5
                   )
lims = unique(lims)
lims = split(lims, lims$Panel) # make a list
lims = lapply(lims, function(l) {
  scale_y_continuous(limits = c(l$ymin, l$ymax),
                     n.breaks = l$n)
  }
)


title = 'Desiree'
data.SE = data.long.SE[data.long.SE$Genotype == title, ]
myggplot(data.SE, title, cnt, lncol = 2, lims, subdir = 'metabolomics')

4.4 correlation

mytitle = 'Desiree'
df = omics[omics$Genotype == mytitle, ]
df = as.matrix(df[, nx:ncol(df)])
df = apply(df, 2, as.numeric)
cs = colSums(is.na(df))/nrow(df)
if (any(cs == 1)) df = df[, -which(cs == 1)]

rownames(df) = omics$SampleID[omics$Genotype == mytitle]

mycond = levels(omics$Condition)


# Open the PDF device
pdf(file.path('..', 'reports', 'metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")), 
    width = 13, height = 13)

# Create the plot
my.pairs.panels(df, mytitle, mycond, cccex = 1)

# Close the device to save the plot
dev.off()
## png 
##   2
# Now plot it to display in the document
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'metabolomics')





mycond = 'C'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'metabolomics')


mycond = 'H'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])



pdf(
         file.path('..', 'reports', 'metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'metabolomics')


mycond = 'D'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'metabolomics')


mycond = 'W'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'metabolomics')

    
mycond = 'HD'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'metabolomics')
genotype = 'Desiree'
De.log2FC = NULL
De.p = NULL

for (i in 1:nrow(comp)) {
  
  # print(i)
  # print( as.character(comp[i, ]))
  
  dtX = data.long[data.long$mygroup %in% comp[i, c(1,3)] & data.long$Genotype %in% genotype, ]

  stress = dtX[dtX$mygroup == as.character(comp[i, 1]), ]
  Ctrl = dtX[dtX$mygroup == as.character(comp[i, 3]), ]
  
  StressC = gsub('_.*', '', comp[i, 1])
  CtrlC = gsub('_.*', '', comp[i, 3])
  
  log2FC = my.logFC(Ctrl = Ctrl, stress = stress, StressC = StressC, CtrlC = CtrlC, genotype = genotype)
  p = my.customised.t.test(Ctrl = Ctrl, Stress = stress, StressC = StressC, CtrlC = CtrlC, genotype = genotype)
  
  De.log2FC = rbind(De.log2FC, log2FC)
  De.p = rbind(De.p, p)
  
  
}


Desiree.stat = merge(De.log2FC, De.p,
                     by.y = c("Genotype", "molecule", "group2", "group1"),
                     by.x = c("Genotype", "molecule", "stressGroup", "controlGroup"),
                     all = TRUE)
Desiree.stat$stressGroup = paste(gsub('_.*', '', Desiree.stat$stressGroup),
                                  stringr::str_pad(gsub('.*_', '', Desiree.stat$stressGroup), 2, pad = "0"),
                                  sep = '_')
logFC = tidyr::spread(Desiree.stat[, c(2:3, 7)], stressGroup, log2FC)
p = tidyr::spread(Desiree.stat[, c(2:3, 13)], stressGroup, p)
colnames(logFC)[2:ncol(logFC)] = paste(colnames(logFC)[2:ncol(logFC)], 'logFC', sep = ' | ')
colnames(p)[2:ncol(p)] = paste(colnames(p)[2:ncol(p)], 'p', sep = ' | ')
stat = merge(logFC, p, by = 'molecule')
ind = match(sort(colnames(stat)[2:ncol(stat)]), colnames(stat)[2:ncol(stat)])
stat = stat[, c(1, ind+1)]
stat[2:ncol(stat)] = format(round(as.data.frame(stat[2:ncol(stat)]), 3), nsmall = 3)
openxlsx::write.xlsx(stat, '../reports/Desiree.logFC-metabolomics_leaves.xlsx', asTable = TRUE)

4.5 heatmap

comp = as.data.frame(comp)
forHeatmap = comp[intersect(grep('C', comp[, 1], invert = TRUE), grep('C', comp[, 3])), ]
forHeatmap$comparison = paste(forHeatmap[, 1], forHeatmap[, 3], sep = '-')

df = Desiree.stat[Desiree.stat$stressGroup %in% forHeatmap[, 1] & Desiree.stat$controlGroup %in% forHeatmap[, 3], ]

my.pheatmap.logFC(df, forHeatmap,
                  main = 'Desiree-metabolomics',
                  cluster_rows = FALSE, cluster_cols = FALSE,
                  subdir = 'metabolomics',
                  width = 6, height = 6, cnt = cnt)
metabolomics.leaves = stat[, grep('molecule|log', colnames(stat))]
# rm()
gc()
gc()

5 hormonomics

dir.create(file.path('..', 'reports', 'hormonomics'))

cnt = 2
fpi = file.path('..', 'input')


fn = "data_targeted.xlsx"
omics = openxlsx::read.xlsx(file.path(fpi, fn), sheet = 'hormonomics-leaves')
data.table::setDF(omics)

Genotype = 'Desiree'
omics = cbind(Genotype, omics)


PostInductionDay = as.numeric(stringr::str_split_i(gsub('S|P', '', omics$SampleID), '_', i = 2))
omics = cbind(PostInductionDay, omics)

omics = omics[omics$PostInductionDay %in% 1:14, ]

Condition = gsub('_.*', '', omics$SampleID)
omics = cbind(Condition, omics)
omics = omics[-which(omics$Condition == 'W' & omics$PostInductionDay %in% c(8,14)), ]

omics$Condition = factor(omics$Condition, levels = c('C', 'H', 'D', 'HD', 'W'))

table(omics$Condition, omics$PostInductionDay)
##     
##      1 7 8 14
##   C  4 4 4  4
##   H  4 4 4  4
##   D  0 0 4  4
##   HD 0 0 4  4
##   W  4 4 0  0
mygroup = paste(omics$Condition, stringr::str_pad(omics$PostInductionDay, 1, pad = "0"), sep = '_')
# table(mygroup %in% comp$V1)
# table(mygroup %in% comp$V3)

omics = cbind(mygroup, omics)

5.1 palette

my.ggplot.palette = c('#7F7F7F', '#C00000', '#FFC000', '#ED7D31',  
                      # '#70AD47', 
                      '#4472C4')

extend.palette = cbind(my.ggplot.palette, c('C', 'H', 'D', 'HD',
                                            #'HDW', 
                                            'W'))
colnames(extend.palette) = c('colour', 'group')
extend.palette = merge(extend.palette, omics[, 1:4], by.x = 'group', by.y = 'Condition', all = TRUE)

5.2 density

# column in which measurements start
nx = 6

data = omics
# colnames(data)[nx:ncol(data)]

data[, nx:ncol(data)] = apply(data[, nx:ncol(data)] , 2, as.numeric)

data.long = tidyr::gather(data, molecule, measurement, colnames(data)[nx]:colnames(data)[ncol(data)], factor_key=FALSE)


data.long$measurement = as.numeric(data.long$measurement)

data.long$molecule = factor(data.long$molecule, levels = colnames(data)[nx:ncol(data)])

data.long$PostInductionDay = ordered(data.long$PostInductionDay, levels = sort(unique(data.long$PostInductionDay )))


data.long = data.long[!is.na(data.long$measurement), ]


# X11()
ggplot(data.long[data.long$Genotype == 'Desiree', ], aes(log2(measurement), fill = Condition, colour = Condition)) +
  geom_density(alpha = 0.1) + 
  ggtitle('Desiree')  +
  theme_classic() +
  scale_colour_manual(name = 'Legend',
                      values = my.ggplot.palette) +
  scale_fill_manual(name = 'Legend',
                    values = my.ggplot.palette) 

# X11()
ggplot(data.long[data.long$Genotype == 'Desiree', ], aes(log2(measurement), fill = mygroup , colour = mygroup )) +
  geom_density(alpha = 0.1) + 
  ggtitle('Desiree')  +
  theme_classic() +
  scale_colour_manual(name = 'Legend',
                      values = unique(extend.palette[, 2:3])[, 1]) +
  scale_fill_manual(name = 'Legend',
                    values = unique(extend.palette[, 2:3])[, 1]
                    ) 

data.long.SE = summarySE(data = data.long, 
                         measurevar="measurement", 
                         groupvars=c("molecule", "Genotype", "Condition", "PostInductionDay", "mygroup"), 
                         na.rm = TRUE)

5.3 plot measurements

cd = sort(unique(data.long.SE$PostInductionDay))

mypallette = my.ggplot.palette


lims = data.long.SE %>% 
  dplyr::group_by(molecule) %>% 
  dplyr::summarise(Panel = molecule, 
                   ymin=min(mean-se),ymax=max(mean+se),
                   n = 5
                   )
lims = unique(lims)
lims = split(lims, lims$Panel) # make a list
lims = lapply(lims, function(l) {
  scale_y_continuous(limits = c(l$ymin, l$ymax),
                     n.breaks = l$n)
  }
)


title = 'Desiree'
data.SE = data.long.SE[data.long.SE$Genotype == title, ]
myggplot(data.SE, title, cnt, lncol = 2, lims, subdir = 'hormonomics')

omics  = omics[, -grep('neoPA', colnames(omics))]
data.long.SE  = data.long.SE[data.long.SE$molecule != 'neoPA', ]
data.long  = data.long[data.long$molecule != 'neoPA', ]

5.4 correlation

mytitle = 'Desiree'
df = omics[omics$Genotype == mytitle, ]
df = as.matrix(df[, nx:ncol(df)])
df = apply(df, 2, as.numeric)
cs = colSums(is.na(df))/nrow(df)
if (any(cs == 1)) df = df[, -which(cs == 1)]

rownames(df) = omics$SampleID[omics$Genotype == mytitle]

mycond = levels(omics$Condition)


pdf(
         file.path('..', 'reports', 'hormonomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## png 
##   2
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'hormonomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'hormonomics')





mycond = 'C'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'hormonomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'hormonomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'hormonomics')


mycond = 'H'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'hormonomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'hormonomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'hormonomics')


mycond = 'D'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'hormonomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'hormonomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'hormonomics')


mycond = 'W'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'hormonomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'hormonomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'hormonomics')

    
mycond = 'HD'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'hormonomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'hormonomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'hormonomics')
genotype = 'Desiree'
De.log2FC = NULL
De.p = NULL

for (i in 1:nrow(comp)) {
  
  # print(i)
  # print( as.character(comp[i, ]))
  
  dtX = data.long[data.long$mygroup %in% comp[i, c(1,3)] & data.long$Genotype %in% genotype, ]

  stress = dtX[dtX$mygroup == as.character(comp[i, 1]), ]
  Ctrl = dtX[dtX$mygroup == as.character(comp[i, 3]), ]
  
  StressC = gsub('_.*', '', comp[i, 1])
  CtrlC = gsub('_.*', '', comp[i, 3])
  
  log2FC = my.logFC(Ctrl = Ctrl, stress = stress, StressC = StressC, CtrlC = CtrlC, genotype = genotype)
  p = my.customised.t.test(Ctrl = Ctrl, Stress = stress, StressC = StressC, CtrlC = CtrlC, genotype = genotype)
  
  De.log2FC = rbind(De.log2FC, log2FC)
  De.p = rbind(De.p, p)
  
  
}


Desiree.stat = merge(De.log2FC, De.p,
                     by.y = c("Genotype", "molecule", "group2", "group1"),
                     by.x = c("Genotype", "molecule", "stressGroup", "controlGroup"),
                     all = TRUE)
Desiree.stat$stressGroup = paste(gsub('_.*', '', Desiree.stat$stressGroup),
                                  stringr::str_pad(gsub('.*_', '', Desiree.stat$stressGroup), 2, pad = "0"),
                                  sep = '_')
logFC = tidyr::spread(Desiree.stat[, c(2:3, 7)], stressGroup, log2FC)
p = tidyr::spread(Desiree.stat[, c(2:3, 13)], stressGroup, p)
colnames(logFC)[2:ncol(logFC)] = paste(colnames(logFC)[2:ncol(logFC)], 'logFC', sep = ' | ')
colnames(p)[2:ncol(p)] = paste(colnames(p)[2:ncol(p)], 'p', sep = ' | ')
stat = merge(logFC, p, by = 'molecule')
ind = match(sort(colnames(stat)[2:ncol(stat)]), colnames(stat)[2:ncol(stat)])
stat = stat[, c(1, ind+1)]
stat[2:ncol(stat)] = format(round(as.data.frame(stat[2:ncol(stat)]), 3), nsmall = 3)
openxlsx::write.xlsx(stat, '../reports/Desiree.logFC-hormonomics_leaves.xlsx', asTable = TRUE)

5.5 heatmap

comp = as.data.frame(comp)
forHeatmap = comp[intersect(grep('C', comp[, 1], invert = TRUE), grep('C', comp[, 3])), ]
forHeatmap$comparison = paste(forHeatmap[, 1], forHeatmap[, 3], sep = '-')

df = Desiree.stat[Desiree.stat$stressGroup %in% forHeatmap[, 1] & Desiree.stat$controlGroup %in% forHeatmap[, 3], ]

my.pheatmap.logFC(df, forHeatmap,
                  main = 'Desiree-hormonomics',
                  cluster_rows = FALSE, cluster_cols = FALSE,
                  subdir = 'hormonomics',
                  width = 6, height = 6, cnt = cnt)
hormonomics.leaves = stat[, grep('molecule|log', colnames(stat))]

# rm()
gc()
gc()

6 transcriptomics

dir.create(file.path('..', 'reports', 'transcriptomics'))

cnt = 3
fpi = file.path('..', 'input')


fn = "data_targeted.xlsx"
omics = openxlsx::read.xlsx(file.path(fpi, fn), sheet = 'transcriptomics-leaves')
data.table::setDF(omics)

Genotype = 'Desiree'
omics = cbind(Genotype, omics)


PostInductionDay = as.numeric(stringr::str_split_i(gsub('S|P', '', omics$SampleID), '_', i = 2))
omics = cbind(PostInductionDay, omics)

omics = omics[omics$PostInductionDay %in% 1:14, ]

Condition = gsub('_.*', '', omics$SampleID)
omics = cbind(Condition, omics)
omics = omics[-which(omics$Condition == 'W' & omics$PostInductionDay %in% c(8,14)), ]

omics$Condition = factor(omics$Condition, levels = c('C', 'H', 'D', 'HD', 'W'))

table(omics$Condition, omics$PostInductionDay)
##     
##      1 7 8 14
##   C  4 4 4  4
##   H  4 4 4  4
##   D  0 0 4  4
##   HD 0 0 4  4
##   W  4 4 0  0
mygroup = paste(omics$Condition, stringr::str_pad(omics$PostInductionDay, 1, pad = "0"), sep = '_')
# table(mygroup %in% comp$V1)
# table(mygroup %in% comp$V3)

omics = cbind(mygroup, omics)

6.1 palette

my.ggplot.palette = c('#7F7F7F', '#C00000', '#FFC000', '#ED7D31',  
                      # '#70AD47', 
                      '#4472C4')

extend.palette = cbind(my.ggplot.palette, c('C', 'H', 'D', 'HD',
                                            #'HDW', 
                                            'W'))
colnames(extend.palette) = c('colour', 'group')
extend.palette = merge(extend.palette, omics[, 1:4], by.x = 'group', by.y = 'Condition', all = TRUE)

6.2 density

# column in which measurements start
nx = 6

data = omics
# colnames(data)[nx:ncol(data)]

data[, nx:ncol(data)] = apply(data[, nx:ncol(data)] , 2, as.numeric)

data.long = tidyr::gather(data, molecule, measurement, colnames(data)[nx]:colnames(data)[ncol(data)], factor_key=FALSE)


data.long$measurement = as.numeric(data.long$measurement)

data.long$molecule = factor(data.long$molecule, levels = colnames(data)[nx:ncol(data)])

data.long$PostInductionDay = ordered(data.long$PostInductionDay, levels = sort(unique(data.long$PostInductionDay )))


data.long = data.long[!is.na(data.long$measurement), ]


# X11()
ggplot(data.long[data.long$Genotype == 'Desiree', ], aes(log2(measurement), fill = Condition, colour = Condition)) +
  geom_density(alpha = 0.1) + 
  ggtitle('Desiree')  +
  theme_classic() +
  scale_colour_manual(name = 'Legend',
                      values = my.ggplot.palette) +
  scale_fill_manual(name = 'Legend',
                    values = my.ggplot.palette) 

# X11()
ggplot(data.long[data.long$Genotype == 'Desiree', ], aes(log2(measurement), fill = mygroup , colour = mygroup )) +
  geom_density(alpha = 0.1) + 
  ggtitle('Desiree')  +
  theme_classic() +
  scale_colour_manual(name = 'Legend',
                      values = unique(extend.palette[, 2:3])[, 1]) +
  scale_fill_manual(name = 'Legend',
                    values = unique(extend.palette[, 2:3])[, 1]
                    ) 

data.long.SE = summarySE(data = data.long, 
                         measurevar="measurement", 
                         groupvars=c("molecule", "Genotype", "Condition", "PostInductionDay", "mygroup"), 
                         na.rm = TRUE)

6.3 plot measurements

cd = sort(unique(data.long.SE$PostInductionDay))

mypallette = my.ggplot.palette


lims = data.long.SE %>% 
  dplyr::group_by(molecule) %>% 
  dplyr::summarise(Panel = molecule, 
                   ymin=min(mean-se),ymax=max(mean+se),
                   n = 5
                   )
lims = unique(lims)
lims = split(lims, lims$Panel) # make a list
lims = lapply(lims, function(l) {
  scale_y_continuous(limits = c(l$ymin, l$ymax),
                     n.breaks = l$n)
  }
)


title = 'Desiree'
data.SE = data.long.SE[data.long.SE$Genotype == title, ]
myggplot(data.SE, title, cnt, lncol = 2, lims, subdir = 'transcriptomics')

6.4 correlation

mytitle = 'Desiree'
df = omics[omics$Genotype == mytitle, ]
df = as.matrix(df[, nx:ncol(df)])
df = apply(df, 2, as.numeric)
cs = colSums(is.na(df))/nrow(df)
if (any(cs == 1)) df = df[, -which(cs == 1)]

rownames(df) = omics$SampleID[omics$Genotype == mytitle]

mycond = levels(omics$Condition)


pdf(
         file.path('..', 'reports', 'transcriptomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## png 
##   2
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'transcriptomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'transcriptomics')





mycond = 'C'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'transcriptomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'transcriptomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'transcriptomics')


mycond = 'H'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'transcriptomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'transcriptomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'transcriptomics')


mycond = 'D'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'transcriptomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'transcriptomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'transcriptomics')


mycond = 'W'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'transcriptomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'transcriptomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'transcriptomics')

    
mycond = 'HD'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'transcriptomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'transcriptomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'transcriptomics')
genotype = 'Desiree'
De.log2FC = NULL
De.p = NULL

for (i in 1:nrow(comp)) {
  
  # print(i)
  # print( as.character(comp[i, ]))
  
  dtX = data.long[data.long$mygroup %in% comp[i, c(1,3)] & data.long$Genotype %in% genotype, ]

  stress = dtX[dtX$mygroup == as.character(comp[i, 1]), ]
  Ctrl = dtX[dtX$mygroup == as.character(comp[i, 3]), ]
  
  StressC = gsub('_.*', '', comp[i, 1])
  CtrlC = gsub('_.*', '', comp[i, 3])
  
  log2FC = my.logFC(Ctrl = Ctrl, stress = stress, StressC = StressC, CtrlC = CtrlC, genotype = genotype)
  p = my.customised.t.test(Ctrl = Ctrl, Stress = stress, StressC = StressC, CtrlC = CtrlC, genotype = genotype)
  
  De.log2FC = rbind(De.log2FC, log2FC)
  De.p = rbind(De.p, p)
  
  
}


Desiree.stat = merge(De.log2FC, De.p,
                     by.y = c("Genotype", "molecule", "group2", "group1"),
                     by.x = c("Genotype", "molecule", "stressGroup", "controlGroup"),
                     all = TRUE)
Desiree.stat$stressGroup = paste(gsub('_.*', '', Desiree.stat$stressGroup),
                                  stringr::str_pad(gsub('.*_', '', Desiree.stat$stressGroup), 2, pad = "0"),
                                  sep = '_')
logFC = tidyr::spread(Desiree.stat[, c(2:3, 7)], stressGroup, log2FC)
p = tidyr::spread(Desiree.stat[, c(2:3, 13)], stressGroup, p)
colnames(logFC)[2:ncol(logFC)] = paste(colnames(logFC)[2:ncol(logFC)], 'logFC', sep = ' | ')
colnames(p)[2:ncol(p)] = paste(colnames(p)[2:ncol(p)], 'p', sep = ' | ')
stat = merge(logFC, p, by = 'molecule')
ind = match(sort(colnames(stat)[2:ncol(stat)]), colnames(stat)[2:ncol(stat)])
stat = stat[, c(1, ind+1)]
stat[2:ncol(stat)] = format(round(as.data.frame(stat[2:ncol(stat)]), 3), nsmall = 3)
openxlsx::write.xlsx(stat, '../reports/Desiree.logFC-transcriptomics_leaves.xlsx', asTable = TRUE)

6.5 heatmap

comp = as.data.frame(comp)
forHeatmap = comp[intersect(grep('C', comp[, 1], invert = TRUE), grep('C', comp[, 3])), ]
forHeatmap$comparison = paste(forHeatmap[, 1], forHeatmap[, 3], sep = '-')

df = Desiree.stat[Desiree.stat$stressGroup %in% forHeatmap[, 1] & Desiree.stat$controlGroup %in% forHeatmap[, 3], ]

my.pheatmap.logFC(df, forHeatmap,
                  main = 'Desiree-transcriptomics',
                  cluster_rows = FALSE, cluster_cols = FALSE,
                  subdir = 'transcriptomics',
                  width = 6, height = 6, cnt = cnt)
transcriptomics.leaves = stat[, grep('molecule|log', colnames(stat))]

# rm()
gc()
gc()

7 tubers metabolomics

dir.create(file.path('..', 'reports', 'tubers-metabolomics'))

cnt = 4
fpi = file.path('..', 'input')


fn = "data_targeted.xlsx"
omics = openxlsx::read.xlsx(file.path(fpi, fn), sheet = 'metabolomics-tubers')
data.table::setDF(omics)

Genotype = 'Desiree'
omics = cbind(Genotype, omics)


PostInductionDay = as.numeric(stringr::str_split_i(gsub('S|P', '', omics$SampleID), '_', i = 2))
omics = cbind(PostInductionDay, omics)


Condition = gsub('_.*', '', omics$SampleID)
omics = cbind(Condition, omics)

omics$Condition = factor(omics$Condition, levels = c('C', 'H', 'D', 'HD', 'W'))

table(omics$Condition, omics$PostInductionDay)
##     
##      28
##   C  12
##   H  12
##   D  12
##   HD 12
##   W   8
mygroup = paste(omics$Condition, stringr::str_pad(omics$PostInductionDay, 1, pad = "0"), sep = '_')
# table(mygroup %in% comp$V1)
# table(mygroup %in% comp$V3)

omics = cbind(mygroup, omics)

7.1 palette

my.ggplot.palette = c('#7F7F7F', '#C00000', '#FFC000', '#ED7D31',  
                      # '#70AD47', 
                      '#4472C4')

extend.palette = cbind(my.ggplot.palette, c('C', 'H', 'D', 'HD',
                                            #'HDW', 
                                            'W'))
colnames(extend.palette) = c('colour', 'group')
extend.palette = merge(extend.palette, omics[, 1:4], by.x = 'group', by.y = 'Condition', all = TRUE)

7.2 density

# column in which measurements start
nx = 9

data = omics
# colnames(data)[nx:ncol(data)]

data[, nx:ncol(data)] = apply(data[, nx:ncol(data)] , 2, as.numeric)

data.long = tidyr::gather(data, molecule, measurement, colnames(data)[nx]:colnames(data)[ncol(data)], factor_key=FALSE)


data.long$measurement = as.numeric(data.long$measurement)

data.long$molecule = factor(data.long$molecule, levels = colnames(data)[nx:ncol(data)])

data.long$PostInductionDay = ordered(data.long$PostInductionDay, levels = sort(unique(data.long$PostInductionDay )))


data.long = data.long[!is.na(data.long$measurement), ]


# X11()
ggplot(data.long[data.long$Genotype == 'Desiree', ], aes(log2(measurement), fill = Condition, colour = Condition)) +
  geom_density(alpha = 0.1) + 
  ggtitle('Desiree')  +
  theme_classic() +
  scale_colour_manual(name = 'Legend',
                      values = my.ggplot.palette) +
  scale_fill_manual(name = 'Legend',
                    values = my.ggplot.palette) 

# X11()
ggplot(data.long[data.long$Genotype == 'Desiree', ], aes(log2(measurement), fill = mygroup , colour = mygroup )) +
  geom_density(alpha = 0.1) + 
  ggtitle('Desiree')  +
  theme_classic() +
  scale_colour_manual(name = 'Legend',
                      values = unique(extend.palette[, 2:3])[, 1]) +
  scale_fill_manual(name = 'Legend',
                    values = unique(extend.palette[, 2:3])[, 1]
                    ) 

data.long.SE = summarySE(data = data.long, 
                         measurevar="measurement", 
                         groupvars=c("molecule", "Genotype", "Condition", "PostInductionDay", "mygroup"), 
                         na.rm = TRUE)

7.3 plot measurements

cd = sort(unique(data.long.SE$PostInductionDay))

mypallette = my.ggplot.palette




title = 'Tubers'
lncol = 2
subdir = 'tubers-metabolomics'

# X11();
myplotY = ggplot(data.long, aes(x = Condition , 
                                y = measurement, 
                                group = Condition,
                                shape = Condition,
                                fill = Condition,
                                colour = Condition)) + 
  geom_boxplot() + 
      facet_wrap(~ molecule, ncol = ncol,
             strip.position = "top",
             scales = 'free_y',
             drop = FALSE
             ) +
      ggtitle(title)  +
    theme_classic() +
      scale_colour_manual(name = 'Legend',
                          values = mypallette) +
      scale_fill_manual(name = 'Legend',
                        values = mypallette
                        ) + 
        labs(x = '', 
             y = 'Relative molecule abundance [+/- SE]') + 
        scale_y_continuous(# trans = log2_trans(), 
                           labels = scales::comma
                           ) +

      theme(axis.text.x = element_text (angle = 0, # 45, 
                                        vjust = 1, 
                                        hjust = 1, 
                                        size = 12),
            axis.text.y = element_text (angle = 0, 
                                        vjust = 1, 
                                        hjust = 1, 
                                        size = 12),
            axis.title.x = element_text (angle = 0, # 45, 
                                        #vjust = 1, 
                                        #hjust = 1, 
                                        size = 12,
                                        face = "bold"),
            axis.title.y = element_text (angle = 90, 
                                        # vjust = 1, 
                                        # hjust = 1, 
                                        size = 12,
                                        face = "bold"),
            axis.text = element_text( angle = 90 ),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            panel.border = element_blank(),
            panel.background = element_blank(),
            legend.position = "bottom",
            legend.background = element_rect(fill=NULL, 
                                             size=0.1, linetype="solid"),
            legend.title = element_text(color = NULL,
                                        size = 12,
                                        face = "bold"),
            legend.text = element_text( size = 12),
            strip.text = element_text(size = 12)) + 
  theme(
        legend.box = "horizontal",
        legend.key.width = unit(2, 'cm'),
        legend.position = c(1, 0),
        legend.justification = c(1, 0)) +
  guides(colour = guide_legend(ncol = lncol)) +
  guides(linetype = guide_legend(ncol = lncol)) +
  guides(shape = "none") +
  theme(panel.spacing.y = unit(1.5, "lines")) +
  labs(color  = "Legend", linetype = "Legend", shape = "Legend")

plot(myplotY)  

# dev.off()
  
  
  fng = paste0(stringr::str_pad(cnt, 3, pad = "0"), '_', title, "_L.pdf")
  ggsave(file.path('..', 'reports', subdir, fng), 
         width = 16, height = 12, units = "in")
  fng = paste0(stringr::str_pad(cnt, 3, pad = "0"), '_', title, "_L.svg")
  ggsave(file.path('..', 'reports', subdir, fng),
         width = 16, height = 12, units = "in")

7.4 correlation

mytitle = 'Desiree'
df = omics[omics$Genotype == mytitle, ]
df = as.matrix(df[, nx:ncol(df)])
df = apply(df, 2, as.numeric)
cs = colSums(is.na(df))/nrow(df)
if (any(cs == 1)) df = df[, -which(cs == 1)]

rownames(df) = omics$SampleID[omics$Genotype == mytitle]

mycond = levels(omics$Condition)


pdf(
         file.path('..', 'reports', 'tubers-metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## png 
##   2
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'tubers-metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'tubers-metabolomics')





mycond = 'C'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'tubers-metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'tubers-metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'tubers-metabolomics')


mycond = 'H'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'tubers-metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'tubers-metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'tubers-metabolomics')


mycond = 'D'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'tubers-metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'tubers-metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'tubers-metabolomics')


mycond = 'W'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'tubers-metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'tubers-metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'tubers-metabolomics')

    
mycond = 'HD'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])


pdf(
         file.path('..', 'reports', 'tubers-metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.off()
## pdf 
##   3
my.pairs.panels(df, mytitle, mycond, cccex = 1)

my.heatmaply(df, mytitle, mycond, subdir = 'tubers-metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'tubers-metabolomics')
genotype = 'Desiree'
De.log2FC = NULL
De.p = NULL

comp = rbind(cbind('H_28',  '-', 'C_28'),
             cbind('D_28',  '-', 'C_28'),
             cbind('HD_28',  '-', 'C_28'),
             cbind('W_28',  '-', 'C_28'))

for (i in 1:nrow(comp)) {
  
  # print(i)
  # print( as.character(comp[i, ]))
  
  dtX = data.long[data.long$mygroup %in% comp[i, c(1,3)] & data.long$Genotype %in% genotype, ]

  stress = dtX[dtX$mygroup == as.character(comp[i, 1]), ]
  Ctrl = dtX[dtX$mygroup == as.character(comp[i, 3]), ]
  
  StressC = gsub('_.*', '', comp[i, 1])
  CtrlC = gsub('_.*', '', comp[i, 3])
  
  log2FC = my.logFC(Ctrl = Ctrl, stress = stress, StressC = StressC, CtrlC = CtrlC, genotype = genotype)
  p = my.customised.t.test(Ctrl = Ctrl, Stress = stress, StressC = StressC, CtrlC = CtrlC, genotype = genotype)
  
  De.log2FC = rbind(De.log2FC, log2FC)
  De.p = rbind(De.p, p)
  
  
}


Desiree.stat = merge(De.log2FC, De.p,
                     by.y = c("Genotype", "molecule", "group2", "group1"),
                     by.x = c("Genotype", "molecule", "stressGroup", "controlGroup"),
                     all = TRUE)
Desiree.stat$stressGroup = paste(gsub('_.*', '', Desiree.stat$stressGroup),
                                  stringr::str_pad(gsub('.*_', '', Desiree.stat$stressGroup), 2, pad = "0"),
                                  sep = '_')
logFC = tidyr::spread(Desiree.stat[, c(2:3, 7)], stressGroup, log2FC)
p = tidyr::spread(Desiree.stat[, c(2:3, 13)], stressGroup, p)
colnames(logFC)[2:ncol(logFC)] = paste(colnames(logFC)[2:ncol(logFC)], 'logFC', sep = ' | ')
colnames(p)[2:ncol(p)] = paste(colnames(p)[2:ncol(p)], 'p', sep = ' | ')
stat = merge(logFC, p, by = 'molecule')
ind = match(sort(colnames(stat)[2:ncol(stat)]), colnames(stat)[2:ncol(stat)])
stat = stat[, c(1, ind+1)]
stat[2:ncol(stat)] = format(round(as.data.frame(stat[2:ncol(stat)]), 3), nsmall = 3)
openxlsx::write.xlsx(stat, '../reports/Desiree.logFC-metabolomics_tubers.xlsx', asTable = TRUE)

7.5 heatmap

comp = as.data.frame(comp)
forHeatmap = comp[intersect(grep('C', comp[, 1], invert = TRUE), grep('C', comp[, 3])), ]
forHeatmap$comparison = paste(forHeatmap[, 1], forHeatmap[, 3], sep = '-')

df = Desiree.stat[Desiree.stat$stressGroup %in% forHeatmap[, 1] & Desiree.stat$controlGroup %in% forHeatmap[, 3], ]

my.pheatmap.logFC(df, forHeatmap,
                  main = 'Desiree-tubers-metabolomics',
                  cluster_rows = FALSE, cluster_cols = FALSE,
                  subdir = 'tubers-metabolomics',
                  width = 6, height = 6, cnt = cnt)
metabolomics.tubers = stat[, grep('molecule|log', colnames(stat))]

# rm()
gc()
gc()

8 logFC results

targeted = rbind(metabolomics.leaves, hormonomics.leaves, transcriptomics.leaves)
targeted = merge(targeted, metabolomics.tubers, by = 'molecule', all.x = TRUE, all.y = TRUE)

data.table::fwrite(targeted, file.path('..', 'output', 'targeted_logFC.txt'), sep = '\t')

9 Session info

# devtools::session_info()
sessionInfo()
## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_United Kingdom.utf8 
## [2] LC_CTYPE=English_United Kingdom.utf8   
## [3] LC_MONETARY=English_United Kingdom.utf8
## [4] LC_NUMERIC=C                           
## [5] LC_TIME=English_United Kingdom.utf8    
## 
## time zone: Europe/Ljubljana
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ggh4x_0.2.8        htmlwidgets_1.6.4  heatmaply_1.5.0    viridis_0.6.5     
##  [5] viridisLite_0.4.2  plotly_4.10.4      RColorBrewer_1.1-3 corrplot_0.95     
##  [9] rstatix_0.7.2      magrittr_2.0.3     ggplot2_3.5.1      scales_1.3.0      
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6      xfun_0.49         bslib_0.8.0       psych_2.4.6.26   
##  [5] lattice_0.22-6    crosstalk_1.2.1   vctrs_0.6.5       tools_4.4.1      
##  [9] generics_0.1.3    parallel_4.4.1    tibble_3.2.1      ca_0.71.1        
## [13] fansi_1.0.6       pkgconfig_2.0.3   pheatmap_1.0.12   data.table_1.16.2
## [17] assertthat_0.2.1  webshot_0.5.5     lifecycle_1.0.4   farver_2.1.2     
## [21] compiler_4.4.1    stringr_1.5.1     textshaping_0.4.0 mnormt_2.1.1     
## [25] munsell_0.5.1     codetools_0.2-20  carData_3.0-5     seriation_1.5.6  
## [29] htmltools_0.5.8.1 sass_0.4.9        yaml_2.3.10       lazyeval_0.2.2   
## [33] Formula_1.2-5     pillar_1.9.0      car_3.1-3         jquerylib_0.1.4  
## [37] tidyr_1.3.1       cachem_1.1.0      iterators_1.0.14  TSP_1.2-4        
## [41] abind_1.4-8       foreach_1.5.2     nlme_3.1-166      tidyselect_1.2.1 
## [45] zip_2.3.1         digest_0.6.37     stringi_1.8.4     reshape2_1.4.4   
## [49] dplyr_1.1.4       purrr_1.0.2       labeling_0.4.3    fastmap_1.2.0    
## [53] grid_4.4.1        colorspace_2.1-1  cli_3.6.3         utf8_1.2.4       
## [57] broom_1.0.7       withr_3.0.2       backports_1.5.0   registry_0.5-1   
## [61] rmarkdown_2.29    httr_1.4.7        gridExtra_2.3     ragg_1.3.3       
## [65] openxlsx_4.2.7.1  evaluate_1.0.1    knitr_1.49        rlang_1.1.4      
## [69] Rcpp_1.0.13-1     dendextend_1.19.0 glue_1.8.0        svglite_2.1.3    
## [73] rstudioapi_0.17.1 jsonlite_1.8.9    plyr_1.8.9        R6_2.5.1         
## [77] systemfonts_1.1.0